﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace DetectScars.common
{
    internal class OtherUtils
    {

        /// <summary>
        /// 生成三段式曲线：常数-样条-常数
        /// </summary>
        public static (double[] xnew, float[] ynew) GeneratePoints(
            double[] x, double[] y, double a, double b, int newCount)
        {
            if (x == null || y == null || x.Length != y.Length)
                throw new ArgumentException("x 与 y 长度不一致或为空");
            if (x.Length < 2)
                throw new ArgumentException("至少需要 2 个点");
            if (newCount < 2)
                throw new ArgumentException("newCount 至少为 2");

            // 1. 生成等距 xnew
            double[] xnew = new double[newCount];
            for (int i = 0; i < newCount; i++)
                xnew[i] = a + (b - a) * i / (newCount - 1.0);

            float[] ynew = new float[newCount];

            // 2. 构造三段
            double x0 = x[0];
            double xn_1 = x[^2];   // C# 8 语法，等价于 x[x.Length-1]

            // 2.1 第一段：常数 y0
            int split1 = Array.FindLastIndex(xnew, v => v <= x0);
            for (int i = 0; i <= split1; i++)
                ynew[i] = (float)y[0];

            // 2.2 第二段：自然样条
            int split2 = Array.FindIndex(xnew, v => v >= xn_1);
            if (split2 < 0) split2 = newCount;   // 防止极端情况

            var spline = new NaturalCubicSpline(x, y);
            for (int i = split1 + 1; i < split2; i++)
                ynew[i] = (float)spline.Interpolate(xnew[i]);

            // 2.3 第三段：常数 yn-1
            for (int i = split2; i < newCount; i++)
                ynew[i] = (float)y[^1];

            return (xnew, ynew);
        }

    }
}
/// <summary>
/// 自然三次样条（边界二阶导为 0）
/// </summary>
public class NaturalCubicSpline
{
    private readonly double[] _x;
    private readonly double[] _y;
    private readonly double[] _h;
    private readonly double[] _a;
    private readonly double[] _b;
    private readonly double[] _c;

    public NaturalCubicSpline(double[] x, double[] y)
    {
        if (x.Length != y.Length) throw new ArgumentException();
        int n = x.Length;
        _x = (double[])x.Clone();
        _y = (double[])y.Clone();

        _h = new double[n - 1];
        for (int i = 0; i < n - 1; i++)
            _h[i] = _x[i + 1] - _x[i];

        // 三对角方程求解二阶导 c[]
        _c = new double[n];
        double[] diag = new double[n];
        double[] rhs = new double[n];

        // 自然边界：c[0]=c[n-1]=0
        for (int i = 1; i < n - 1; i++)
        {
            diag[i] = 2 * (_h[i - 1] + _h[i]);
            rhs[i] = 3 * ((_y[i + 1] - _y[i]) / _h[i] - (_y[i] - _y[i - 1]) / _h[i - 1]);
        }

        // 追赶法
        double[] alpha = new double[n];
        double[] beta = new double[n];
        for (int i = 1; i < n - 1; i++)
        {
            double denom = diag[i] - _h[i - 1] * alpha[i - 1];
            alpha[i] = _h[i] / denom;
            beta[i] = (rhs[i] - _h[i - 1] * beta[i - 1]) / denom;
        }

        for (int i = n - 2; i >= 1; i--)
            _c[i] = beta[i] - alpha[i] * _c[i + 1];

        // 计算 a, b
        _a = new double[n - 1];
        _b = new double[n - 1];
        for (int i = 0; i < n - 1; i++)
        {
            _a[i] = (_c[i + 1] - _c[i]) / (3 * _h[i]);
            _b[i] = (_y[i + 1] - _y[i]) / _h[i] - _h[i] * (_c[i + 1] + 2 * _c[i]) / 3;
        }
    }

    public double Interpolate(double t)
    {
        int n = _x.Length;
        if (t <= _x[0]) return _y[0];
        if (t >= _x[^1]) return _y[^1];

        int idx = Array.BinarySearch(_x, t);
        if (idx < 0) idx = ~idx - 1;
        idx = Math.Max(0, Math.Min(idx, n - 2));

        double dx = t - _x[idx];
        return _y[idx] + dx * (_b[idx] + dx * (_c[idx] + dx * _a[idx]));
    }
}
// 工具：在数组中找满足条件的第一个/最后一个索引
internal static class ArrayExtensions
{
    public static int FindFirstIndex<T>(this T[] array, Predicate<T> match)
    {
        for (int i = 0; i < array.Length; i++)
            if (match(array[i])) return i;
        return -1;
    }

    public static int FindLastIndex<T>(this T[] array, Predicate<T> match)
    {
        for (int i = array.Length - 1; i >= 0; i--)
            if (match(array[i])) return i;
        return -1;
    }
}